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The phase diagram of a half-filled hard core boson two-leg ladder in a flux is investigated by means 
of numerical simulations based on the Density Matrix Renormalization Group (DMRG) algorithm 
and bosonization. We calculate experimentally accessible observables such as the momentum dis¬ 
tribution, as well as rung current, density wave and bond-order wave correlation functions, allowing 
us to identify the Mott Meissner and Mott Vortex states. We follow the transition from commen¬ 
surate Meissner to incommensurate Vortex state at increasing interchain hopping till the critical 
value [Piraud et al. Phys. Rev. B 91 , 140406 (2015)] above which the Meissner state is stable at 
any flux. For flux close to n, and below the critical hopping, we observe the formation of a second 
incommensuration in the Mott Vortex state that could be detectable in current experiments. 

PACS numbers: 03.75.Lm,05.30.Rt,64.70.Rh,71.10.Pm 


Superconductors in external magnetic field H < H c 1 
exhibit the Meissner-Ochsenfcld effect where surface cur¬ 
rents screen completely the magnetic field in the bulk, re¬ 
sulting in perfect diamagnetism |l| • Type-I superconduc¬ 
tors return to the normal state for H > H c 1 while in type- 
II superconductors, for H c 1 < H < H c 2 a vortex phase is 
formed, in which the magnetic field partially penetrates 
the system along flux lines surrounded by screening cur¬ 
rents. This behavior can be understood in the frame¬ 
work of spontaneous breaking of a global U(l) symmetry 
via the Landau-Ginzburg equation[l|. In a quasi-one di¬ 
mensional system, such symmetry breaking is precluded 
by the Mermin-Wagner-Hohenberg theorem^, 1]. How¬ 
ever, in the case of a bosonic two-leg ladder SSEl, 
an analog of the Meissner phase was predicted to exist 
in the ground state for low flux, while for higher flux 
a Tomonaga-Luttinger liquid (TLL) of vortices was ex¬ 
pected. The quantum phase transition between these 
two states is in the commensurate-incommensurate (C- 
IC) universality class [lil,[3. Other orderings have been 
predicted, such as chiral superfluid order at half a flux 
quantum pe r plaqu ette [a, ,91, 1Q| and a chiral Mott insu¬ 
lating phase 11-3], which is a Mott regime [3] possess¬ 
ing chiral currents as well as a spin-density wave phase. 
DMRG studies of ladders with diagonal interchain hop¬ 
ping are also available 17-3]. While the original pro¬ 
posal was made in the context of Josephson junction lad¬ 
ders, where the quantum effects are spoiled by dissipa¬ 
tion [3]) the advent of ultracold atomic gases offers an¬ 
other realization of strongly interacting one dimensional 
boson systems [ 3 . 13 ] . Moreover, it has been shown the¬ 
oretically 26|, |27| and experimentally 


how artificial 

gauge field could be created in these systems. Recently, 


the transition from Meissner to Vortex phases in non¬ 
interacting bosonic ladders of ultracold atoms has been 
studied experimentally at fixed flux tt/2 per plaquette 
and variable interleg hopping (3 |- 

In this Letter we explore the phase diagram of hard¬ 
core spinless bosons on a two-leg ladder at half-filling 
as a function of flux and interchain hopping by means 
of numerical simulations using DMRG algorithm and 
bosonization. We find, in agreement with [301 ]. that hard¬ 
core constraints cause a significant enlargement of the 
Meissner phase over the vortex one with respect to the 
non-interacting case: above a critical value of the in¬ 
terchain hopping [3 the system remains in the Mott- 
Meissner (MM) state for any flux (see Fig. [TJ. Below 
the critical interchain hopping, both the behavior of the 
momentum distribution and of the rung current, show 
that the transition from Mott-Meissner (MM) to Mott- 
Vortex (MV) state falls in the universality class of the 
C-IC transitionpjJ- For fluxes close to 7r, we observe an¬ 
other incommensuration, whose origin is discussed within 
bosonization. 


We consider 14] a two-component system of hard core 
bosons on two leg ladder, with a flux per plaquette A and 
interchain hopping Cl: 


H x = -tJ2(b]^e lXa b J+ i, ct +H.c.) 

+ °E( 6 i,t 6 hl+ H -C-), (!) 


with foj a (bj a ) bosonic creation (annihilation) operator at 
site j, a = ±1/2 the chain index, and te lXa the hopping 
amplitude along the chain a. This Hamiltonian can be 
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mapped onto a system of spin-1/2 [311 ] bosons with spin- 
orbit coupling in a transverse magnetic field [32] with 
each spinor state corresponding to one leg of the ladder. 
For half-filling, i.e. for one boson per rung, at A = 0 and 
^ 0 the ground state of dT[) is a rung-Mott Insulator 5 j. 
For A > 0, according to the bosonization treatment 31, 
two phases with a charge gap are expected [H, [13, EH, 
EH , the Mott-Meissner (MM) and the Mott-Vortex (MV) 
state. In the MM state, for 0 < A < A c , two currents 
of opposite sign flow along the legsj30|, the interchain 
current 


J r (l) = iQ (b] f bi ti , - 


( 2 ) 


has zero expectation value and exponentially decaying 
correlations, and the screening current, i.e the difference 
between the currents of the two legs 

J s = -itY^ {^e lXa b] a b j+h<7 - a e ~ lX,J b) +l rT bj^ , (3) 
ho¬ 


is a smooth function of the applied flux (increasing lin¬ 
early at small flux). On increasing the flux A > A c (f2), 
the system enters the MV state, there is a sudden 
drop|30| of the screening current J s and simultaneously 
the rung current correlations decay becomes algebraic [30j 
with an incommensurate modulation of wavevector q( A). 
Close to the transition point A c (fl), the wavevector 
q( A) ~ \J A 2 — A^. In the non-interacting case, the 
Hamiltonian Eq. © can be readily diagonalized [1J, 116] 
and Ac (H) = 2 arctan[f2/(4 1)\. The occurrence of the 
MV phase can be seen out also in the total, as well as, in 


the spin resolved momentum distribution 
tem: 


of the sys- 


n(k) = n a (k) 

a 


L—l 


( 4 ) 
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ground-state energy is of order 5 x 10~ 5 at its most. We 
further extrapolate in the limit m —> oo all the quan¬ 
tities calculated to characterize the phase diagram. In 



a) 



FIG. 1: (color online) Left panel: Phase diagram of © as a 
function of flux per plaquette A and f 1/t. The phase bound¬ 
ary between the Meissner and the Vortex phase is shown by 
the black solid line, displaying the persistence of the Meissner 
phase[3Cj above the threshold > Q c for all fluxes A, except 
A = 7 r. For comparison the red-dashed line shows the bound¬ 
ary A/ 0> (f2) for the non-interacting system. In the shaded 
area a second incommensuration appears. In the green area 
(region II) extra peaks at k = ix develop in the FT of the rung 
current correlations and they become the dominant correla¬ 
tions in the blue region (III). The double line (green vs dark 
red) at A = 7r represents the transition to a localized phase. 
In the right panel we show intensity plots of n(fc) versus A 
and k. In panel (a), at fl/t = 1.75 the system is always in the 
MM phase, indeed only a single maximum at k = 0 is visible 
for all A. At A = n, n(k) = 1 indicating the formation of a 
fully localized state (dark red solid line). In panel (b) and (c), 
for SI /t = 1 and 0.25 respectively, the transition from MM to 
the MV with two maxima symmetric around k = 0, is shown. 


In the MM phase n(k) has a single maximum at k = 0, 
whereas in the MV phase it exhibits a pair of max¬ 
ima k = ±q( A)/2[i3]. We have obtained the ground 
state phase-diagram of © by computing various observ¬ 
ables like the momentum distribution and the screen¬ 
ing current J s together with the Fourier Transform (FT) 
C(k) = Yli e~ lkl {J r (l)Jr{ 0)) of the rung current correla¬ 
tion function. 

While performing simulations with both periodic (PBC) 
and open (OBC) boundary conditions, we found the for¬ 
mer to be more suitable for our system, despite the well- 
known computationally more demanding convergence 
properties typical of PBC [35l - l37[ . As such we run sim¬ 
ulations employing PBC for system sizes ranging from 
L = 16 to L = 64, keeping up to m = 1256 states during 
the renormalization procedure. In this way the trunca¬ 
tion error i.e. the weight of the discarded states, is at 
most of order 10 6 , while the maximum error on the 


Fig. m we summarize our findings for the phase diagram 
at half-filling. At variance with the non-interacting case 
where there is a critical A^^H) for all H, in the pres¬ 
ence of the hard-core interaction, for interchain hoppings 
> f2 c , the commensurate-incommensurate transition 
disappears[30j and the MM phase is stable for all fluxes. 
Another effect of the hard-core interaction, as we will 
discuss below, is that in the Vortex phase, at A = 7r and 
A close to 7 r, a commensurate peak appears in C{k ~ 7r), 
along with an incommensuration in the density correla¬ 
tions. At A = 7T, and for H > fl c a fully rung localized 
phase is obtained. Such rung localized ground state was 
discussed in the limit fl t in [30| . 

We have characterized the nature of the Mott-Meissner 
and Mott-Vortex phases by examining C(k ), the stag¬ 
gered boson density wave S (fc) and the symmetric bond- 
order wave Sg OW static structure factors which bring 
information on the spin density and bond-order waves, 
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respectively: 


L — l 


S(k) = — e zfe(j i) sign(fTO- , )(n.,>«;<T'), (5) 


S' 


3 , 1 = 0 


1 L_1 

Bow(k) = y" X/ e ' 

j,i —o 


W-ViSBjSBt); 


( 6 ) 


where Bj = X),, &j+i,A<r + H.c. and SB 3 = B 3 - (£,). 

In Fig. [2] we follow the the MM-MV phase transition 
at small A and fi (see cut one in Fig.l). As predicted 
from bosonization ;32| the vortex phase is signalled by 
the appearance in C(k) of two cusp-likc peaks respec¬ 
tively at k = q( A) and k = 2n — q( A) (see panel a) of 
Fig. [2] whose heights do not scale with the size of the 
system (see Fig. 1 of [32J). In MV phase, the spin re¬ 
solved momentum distribution n a {k) shows a symmetric 
peak centered at k = crq(X ), as predicted by bosoniza¬ 
tion. In this region of parameter space the correlation 
length associated with the Mott gap| is comparable 
with the system size, and the peak takes a cusp-like shape 
as in a Tomonaga-Luttinger liquid (3lj, instead of the 
typical Lorentzian-shape expected for a Mott-insulator. 
Also S(k) shows the expected low momentum behaviour 



FIG. 2: We show FT of correlation functions as from DMRG 
simulation for L = 64 at A = 7r/4 for two different values 
of the Q,/t = 0.0625 and 1, respectively in the Vortex (black 
solid line) and Meissner phase( red solid line). Panel a) shows 
the FT of the rung-current correlation function C(k), panel 
b) the spin correlation functions S(k) and panel c) the charge 
bond-order correlation function Sg Qw . In panel d) the spin 
resolved momentum distribution is shown, with n_ CT (fc) = 
n CT ( — k). 


according to bosonization approach: in the MM phase 
S(k) = S( 0) + ak 2 + o(fc 2 ), with S( 0) > 0, while in the 
MV phase S(k) = Ks ^ + o(fc), with K* = 1 (as expected 
for a hard-core boson system) a signature of a TLL of 


vortices. The transition is also seen in S(k ~ n). In the 
MM phase, S(k ~ n) shows a Lorentzian-shaped peak 
while in the MV phase this peak takes a cusp-like shape. 
A similar change across the MM-MV transition is also 
seen in the correlation function Ssowc{k ~ 7r)(see Sup¬ 
plemental Material HI). This description breaks down 



FIG. 3: We show FT of correlation functions as from DMRG 
simulations for L = 64 at A = 7r at various fl/t. Panel a) 
shows the FT of the rung-current correlation function C(fc), 
panel b) the spin correlation function S(k) and panel c) the 
charge bond-order correlation function Sg 0w . In panel d) 
the spin resolved momentum distribution is shown. Red solid 
curves are for fl/f = 1.75 in the fully localised state, while 
black, blue, green and magenta solid lines are respectively for 
n/t = 1.5,1.25, 1 and 0.5. 


when A is no longer a small quantity as q( A) would be 
comparable to the momentum cutoff. 

At A = 7r the major changes from the conventional C-IC 
transition at small flux are observed. To derive the low 
energy Hamiltonian it becomes necessary to choose the 
gauge with the vector potential along the rungs of the 
ladder, so that the interchain hopping reads: 

H hop . = nJ2(-y b l A— (7) 

3, <* 


After applying bosonization, the hopping Hamiltonian 
can be rewritten in terms of a free boson (f) c describ¬ 
ing the total density fluctuations coupled to SU(2) 1 
Wess-Zumino-Novikov-Witten (WZNW) currents Jr^l 
describing the chain antisymmetric density fluctuations 
by a term oc Q cos \/2<j) c (./)( + J|) ( see [32| for details). 
Such a term can be treated in mean-field theory 2C 


This procedure leads to an effective Hamiltonian with a 
gap A c ~ Q 2 for the total density excitations, while the 
antisymmetric density modes remain gapless and develop 
an incommensuration of wavevectorp(fl) oc fl 2 (see Fig. 2 
in HI}). The presence of this predicted incommensura¬ 
tion is visible in the low momentum behaviour of S(k) 
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and C(k ) (panel a) and panel b) of Fig. [3] that become 
°c p(fi)| + |fc+p(n)|), i.e. constant for |fc| < p(Q) 

and linear in k for |fc| > p(Q). In the S^ ow (k) we ob¬ 
serve a cusp at the same vectors p(fi) (panel c) of Fig. [3]). 
As expected, all these correlation functions also develop 
a peak at k = ir/a. A sign of the incommensurability 
at A = 7 r should be visible also in the momentum distri¬ 
bution n a {k) (see Fig. [3l panel d)). In this case, a cal¬ 
culation based on non-abelian bosonization and operator 
product expansion, would lead to three Lorentzian-like 
peaks centered in 7r/(2a) and n/(2a) ±p(Q)/2. However, 
these peaks cannot be separated if the correlation length 
in real space u c /A c ~ Q~ 2 is shorter than the wavelength 
2n/p(Q) ~ B -2 . In the numerical simulations, at L = 64 
in PBC, (see Fig. [3]) a broad peak is observed for k = j-. 

When A < 7 T (second cut in Fig.l) we can proceed 
analogously to the previous case and choose a gauge such 
that: 

H = + H.c.) 

3 ,° 

( 8 ) 


and define <5A = (A — 7 r), so that the bosonized Hamilto¬ 
nian contains the extra term dA(JJ) — J£). For this case, 
the Fourier transform of the rung current correlation will 
present peaks at k = ^ and k = ^ ± \Jp{Q) 2 + (<5A/a) 2 . 
When <5A is increased, these last two peaks become domi¬ 
nant, and we crossover to the behavior already discussed 
for weak A. At A < 7 r the C(k) (see Fig. [4]) shows, beside 
the peak at k = two peaks symmetric around k = — , 
in real space these last two oscillations exibit an expo¬ 
nential decay for f l/t > 1 and a power law for Q/t < 1 
(region III and II in Fig.l). The situation is reversed for 
the oscillation at k = —. At Q/t = 1 all oscillations, for 
systems with L = 64 in PBC, exhibit power law decay. 
The effect of this incommensuration can also be followed 
in the behaviour at small k of S(k) that instead of being 
a constant value for k < ^/p(Q) 2 + (SX/a) 2 shows a lin¬ 
ear behaviour. In Sg OW (k), for Q/t < 1 two symmetric 
peaks are present at k = ±q( A). We checked that the 
phase is a single component Tomonaga-Luttinger liquid 
by computing the Von Neumann entropy at Q/t = 1.5 for 
A = 0 . 757 T and 0.81257T and obtaining the expected log¬ 
arithmic dependence with system sizej30j. ruling out a a 
Chiral Mott insulator 13), [l5[ for A < 7 r. At general com¬ 
mensurate filling n and flux A = 2irn, the incommensura¬ 


tion generating term becomes — iQe 1 ^ 2 ^ 0 ( J^+ J^)+H.c., 
leading to density wave phases with incommensuration. 
Let us finish by noting that such incommensuration is 
specific of hard core boson systems. With less repulsive 
interactions, the term that gives rise to the vortex lattice 
state would be relevant [ 13 , |15|, |42|, while stronger repul¬ 
sion would make the term stabilizing the checkerboard 
density wave relevant. Adding a nearest neighbor intra- 


C/3 



FIG. 4: We show FT of correlation functions as from DMRG 
simulations for L = 64 at A = 37t/4 at various f l/t. Panel a) 
shows the FT of the rung-current correlation function C(k), 
panel b) the spin correlation function S(k) and panel c) the 
symmetric bond-order correlation function Sbowc ■ In panel 
d) the chain resolved momentum distribution is shown, with 
n- a (k ) = n a {—k). Red solid curves are for Q/t = 1.5 in the 
Meissner phase, while black, blue, green and magenta solid 
lines are respectively for f l/t = 1.25,1 and 0.625. 


chain interaction V to the hard core repulsion, a first 
phase transition at V < 0 will separate the vortex lattice 
from the incommensurate state, and a second transition 
at V > 0 will separate the incommensurate state from 
the density wave state. 

In conclusion, we have studied a two-leg hard core 
boson ladder in an artificial gauge field. In contrast to 
the non-interacting case, the vortex phase is suppressed 
when the interchain hopping exceeds a threshold value, 
as found in [30j]. At flux 7 r per plaquette and f \/t > 1.5 
the ground state becomes a tensor product of singly 
occupied rungs, as was expected^] in the Cl/t —> oo 
limit. For Q/t < 1.5, we have obtained an incommen¬ 
surate insulating state similar to the spin-nematic state 
of frustrated XXZ spin chains (UIHI. In the case of a 
system of weakly coupled ladders, a long range ordered 
phase could form in which density wave or rung current 
would possess a long range commensurate order, but 
exponentially damped incommensurate correlations 
would still be present. The presented results could be 
detectable in current experiments with cold atoms [ 20 | 
and the evidence of a persisting Meissner state could be 
relevant for quantum computing purposes in defining a 
stable flux qbit 43, 441- 
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Quantum Phenomena: from strongly correlated system 
to quantum simulators”. 
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Supplemental Material for “Persisting Meissner state and incommensurate phases of 

hard-core boson ladders in a flux” 

Supplementary material for “Incommensurate phases of a hard-core boson two-leg 

ladder in a flux” 


CORRELATION FUNCTIONS AND OBSERVABLES IN BOSONIZATION 


Bosonized Hamiltonian and observables 


Let us first consider a single leg a in the case of 0 = 0, A = 0. Hard core bosons are mapped on non¬ 
interacting spinless fermions by the Jordan-Wigner transformation [sij . These non-interacting spinless fermions are 
then bosonized [S2j|. The bosonized form of the Hamiltonian Hu reads: 


H w = J2 I k^fdflb) 2 + 


27r . 


u 

K ' 


(SI) 


where [(/> CT (x), W a {x')\ = i8 aa '8{x — x ') and it f x n a = 6 a . in Eq. (IS3l) . u is the velocity of excitations, while K is 
the Tomonaga-Luttinger (TL) parameter. For non-interacting spinless fermions, K = 1 and u = 2ta sin(7m/2), where 
n = (n-f + nj.) is the average number of bosons per site and a is the lattice spacing. At half-filling (i. e. n = 1) u = 21. 
The boson annihilation operators are represented asf 


bi „ = e ie ^ a) 


Aq + ^ A n 


o 2im(0(ja)—7r(n (T )x/a) 


m^O 


(S2) 


In the presence of a vector potential along the legs of the ladder, the lattice Hamiltonian can be brought back to the 
case of A = 0 by the canonical transformation bj }lT = e~A x<7 bj a . The bosonization technique can then be applied to 
the Hamiltonian written in terms of the bj,„ bosons. One finds a Hamiltonian of the form (IS1I) with 6 a , (f> a replacing 
Using Eq. (IS2jl . one obtains 6(x) = 6{x) — A ax/a, and <j>{x) = <f>(x) giving the bosonized Hamiltonian 


*«-£/§ 


uK ( ttHcj + a— 
a 


|(d c <b) 2 


Actually, it is convenient to turn to symmetric (c) and antisymmetric (s) representation, (f> c s = , 9 C S = 

so that the full Hamiltonian reads: 


(S3) 


6 t ±0 l 


H = H c + H, 
dx 
2ir 


H r = 


H x = 


f dx 

J 27T 


c Af c (7rn c ) 2 + +-^{d x 4> c ) 2 

An 


u s K s ( ?rn s + —+ -jridxfisY 


(54) 

(55) 

(56) 


where u c K c = u s K s = uK , u c /K c = u s /K s = u/K , and we have used a = ±1/2. 
The boson annihilation operators become: 


b ia = e ^ Ba+2ae A 


'y ' gi™(%/2<±+2o"\/20s — 2ir(n cr )x/a) 


Therefore, 


(b j<7 bl a ) ~ ( e ^ u A e -^s c (o)^ e ^e s Ua) e -^e s (oy + _ _ 


(S7) 


(S8) 


Using instead the bosonized expression for the density p a = _ d x (f> a + Bi sin(2 <j> a — 2it(n a )x/a) one obtains the 

the total boson density ptot = p-f + Pi as: 


1 \/2 

Ptot{x ) =- -d c (j) c ± 2 B 0 e ln “ sin y/2<f> c cos \pl<.j> s 

a 7r 


(S9) 
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and the antisymmetric density a z = (p^ — p \)/2 as: 


\x) = - ^-=d x cj> s + B 0 e™* cossin V2<j> a , 

7TV 2 


(S10) 


where we have assumed (n a ) = 1/2. Now, if we turn on the interchain hopping f2, we will obtain with the help of 
the following bosonized form: 


H t , 


2n 


dx cos V29 s [A(j + 2 Af cos V8(j> c + 2 A\ cos v / 8/> s 


(Sll) 


where ... stands for less relevant operators. It is important to note that in a renormalization group calculation, 
operators cos y/8(j) c and cos-\/8 (f> s are gen erated in the flow. Since we are in the half-filled case, (n-f + nj.) = 1, 


umklapp scattering cos V8(f> c is present [S4 IS5| and opens of a gap in the c modes. 

For (n-|-) = (nj.), the term cos V8(p a is present, but is marginally irrelevant for repulsive interactions. 

We have in the ground state {<j> c ) = 0. Using Eqs. (ISffll . the staggered and uniform components of the total density 
have exponentially decaying correlations. Meanwhile, the antisymmetric density (IS lOh has a simplified expression for 
distances much larger than u c /A c where one can replace cos \[2<\> c by its average. 

Besides density waves, the system can also present bond ordering. The bond order wave order parameter O bow for 
spin a boson is defined by: 




bj<jbj+ 1,<7 — T{ja) + (-yO° BOW (ja), 


where T is the kinetic energy density, and in bosonization 


0 B ow — Cocos( 2^ ct ) 


We can define the two order parameters: 


O c BO w — 'y ' Obow — 2Co cos y/2(f> c cos \Pl<\) s 

a 

°BOW = Y Sign((7)Og 0ff = 2Cp sin V2cj> c sin V2(j) s 


(512) 

(513) 

(514) 

(515) 


In a Mott phase, O s BOW is always short range ordered, while O bow ~ cos-\/2 (j> 8 . If we consider the real space 
Green’s functions for the bosons, due to the long range order of <p c , the exponentials e I/3Sc of the dual field are short 
range ordered, and the boson Green’s functions decay exponentially. Excitations of the total density are solitons and 
antisolitons of topological charge </> c (+oo) — 4> c {— oo) = ±7 t/\/ 2. Such an excitation corresponds to a change of particle 
number ±1, and for fixed particle number density excitations are formed of soliton/antisoliton pairs. 

For A = 0, the term 2flA^ cos V2d s is a relevant perturbation with scaling dimension 1/(2 K s ), opening a gap 
A s ~ £pK s /(AK s -i) j n antisymmetric sector. We have in the ground state (d a ) = ir/y/2 and as a result all the 
fields e l P't >a are short range ordered. Therefore, in the Mott-Meissner phase, all the density wave and bond order 
wave order parameters are short range ordered. Because of the long range order of 9 S we expect that the correlation 
functions {bj,crbj _ a ) are non-zero. A more precise estimate of the gap (albeit without log corrections) can be obtained 
from the results of [S6 5 , £7]. Using these results, we predict that the soliton mass A s behaves as: 


2 r 


A, = 


(sic s — 2 ) ( 1 4iy,) 


QAqCl 


2 K s 

aK s — 1 


^(^t) \ r(^) 


(sie) 


and the soliton/antisoliton dispersion is E s (k) = yj ( u s k ) 2 + A^. In the case of hard core bosons , we hav e to set 
K s = 1 in Ea. dSlbl) . In such case, besides the solitons and antisolitons, there are two breathers [S8HS10j, a light 
breather of mass A s and a heavy breather of mass y/3A s . The topological charge of the 9 S solitons/antisolitons is 
0 s (+oo) — 9 s (—oo) = ±ttV2. The solitons therefore carry a spin current u s K s . In the case where one is considering 
the gap between the ground state and an excited state of total spin current zero (i. e. containing at least one soliton 
and one antisolito n), th e measured gap will be 2A S . The amplitude Aq can be estimated for hard core bosons in the 
case of half-filling. Sll 


A 2 = 2 l/6 e l/2 A —6 ; 


(S17) 
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where A ~ 1.282427 is Glaisher’s constant. Using (IS 161) with K s = 1 and (IS17I) . we find: 

A u s 2r(l/6) /2 1 / 6 e 1 / 2 j4 _6 r(3/4) 7rfia\ 2/3 
s “ T v^r(2/3) V rwi) u, 7J ’ 


(S18) 


for half-filling. Using u s = 2fa, we finally have A s /t = 3.3896(1 1/f) 2 / 3 . The marginally irrelevant operator cos y/8(f> s 
can give rise to logarithmic corrections to that scaling |S12l |S13I | of the form A s ~ fl 2 / 3 1 lnflj 1 / 6 . 


Commensurate Incommensurate transition 


Neglect ing the m arginally irrelevant term cos as in (HH , the Hamiltonian ([S6]) (IS11I) describes the C-IC 


transition 


sB®. When A exceeds the threshold A c ~ (p.A^a/u s ) 2 1 A 2Ks \ it becomes energetically favorable to 


populate the ground state with a finite density of solitons of the field 0 S to form a Tomonaga-Luttinger liquid (TLL) 
of solitons. The low energy properties of that TLL are described by the effective Hamiltonian: 


H* = 


dx 

2tt 


u* s (X)K* s (X)(nn* s ) 2 + -^(d^) 2 


(S19) 


and we have 7rn s (:r) = 7rn*(:r) — q( A)/ a/2; 0 s (x) = 6*{x) — q( X)x/\/2. Near the transition point A c , we have 
q( A) ~ C\/X — A c . Moreover, as A — > X c + 0, K*( A) goes to a limiting value such that S17 . S18l | the scaling 


dimension of cos-\/2 9 S becomes 1. Since the scaling dimension of cos-\/2 6 S wi th a Hamiltonian of the form (IS19I) is 
l/[2Lf*(A)] we have K*{X —>• A c + 0) = 1/2. Using a fermionization method [S14j |. an explicit form of q{X) can be 
obtained for K s = 1/2. 

The antisymmetric leg current (or screening current) operator 


J* (j) = { ael 


Act it 




— tie 


— iAtri t 


u j-\-l,cr u 3i a J ’ 


(S20) 


is obtained by differentiating the Hamiltonian with respect to the parameter A. One finds: 


J s (x) — 


u s K s 

nV2 


rrl I, 


A 

iV2 


In the Meissner phase, (n s ) = 0 so that: 




(S21) 


(S22) 


In the vortex phase, J s = u s K s (X — sign(A)g(A))/(27r). By fermionization (S14l |. q( A) = yA 2 — A 2 is obtained for 
K s = 1 / 2 . 

In the Mott-Vortex state, we obtain the rung current as: 


j±(x) = UHq sin[\/20* (x) — q(X)x]. 


(S23) 


In both the Meissner and the vortex phase, its expectation value (j±) = 0 vanishes. In the Meissner phase, the 
conversion current correlation function C(x) = (j±{x)j±(0)) decays exponentially whereas in the vortex phase: 


C(x) = i(IL4 2 ) 2 


For K* = 1/2, the Fourier transform: 


C(k) = 


a(^A 2 ) 2 ^_ lk _ q(x) \ a 


cos[< 7 (A)x] 


Q -|fe+(?(A)|a 


), 


(S24) 


(S25) 


has two cusps at k = ±</(A). Peaks divergent with system size appear in C(k) for K* > 1. 
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We can also obtain the spin-spin and bond order correlation functions. In the Mott-Vortex phase, the fields e 
have quasi-long range order. As a result, we find that: 

2 \ K ! 2 

a x ‘ 


(Obow(. x )Obow(®)) — Cq ((cos V2c/) C )) 2 


x 2 + a 2 


(<t z (x)<j z (0)) = ^ 7 ^^ + (-)fB 0 2 ((cosv / 2^}) 2 


47 t 2 (x 2 + a 2 ) 2 

If we turn to the Fourier transforms,we find that for k ~ 0 


K*J2 


S(k) = ^^e- |fc|a , 


(526) 

(527) 

(528) 


by using the integral ^ e +a -i = fe For k ~ ^ and AT* < 1, we find that the correlation functions S(k) and 
Ssowcik) are divergent as: 


S(k) 


'bow 


(k) 


k- 


k:-i 


(S29) 


and as a result a divergence going as \ka\~ 1 ^ 2 is expected at the transition, while far from the transition K* ~ 1, 
giving only a weak power law or logarithmic divergence. For K* > 1, both S s (k) and for S£ ow (k) remain finite in the 
vicinity of tv/ a. 


o; 

o 

8 

g 

g 

O 


CvJ 

g 

« 



FIG. SI: Size dependence of S(k = it) (red dots) and Sgow(k = tt) (blue dots), for system sizes L = 16,32,48,64 and 96. 
Data are for f2 = 0.25 and A = 7r/2, well inside the vortex phase region. Both quantities don’t show a visible size dependence 
in agreement with (IS29I) for K* = 1. Solid black dots are C(k = q(\))a with a = ff 2 /16 to be on the same scale on the graph. 
The size-scaling is compatible with a logarithmic dependence (red dotted line) as expected far from the transition region 


If we turn to the momentum distribution, we will find: 

(bj a bJ ,) = l 5 (T(T /e _iCT ' 3 ( A ) a ( ;,_i ) ^ e i6c ^ a ^'^ 2 e~ iec( - la ^^ 2 )((e ie ^ : ’ a ^^e~ ies( - la ^^ 2 ) 
... . ( e i6 c (ja)/V2 -ie c (la)/V2\ 


= ^ 


(1 + (x/a) 2 ) 


2 a/(8J<rj) 


(530) 

(531) 


because of the exponential decay of the charge correlator, we expect that the momentum distribution of spin er particles 
will be centered around k = — crq(X). The total momentum distribution will thus have two peaks for k = ±g(A)/2. 


INCOMMENSURATION FOR A ~ tt 


For A close to 7r, the form (IS6I) for the Hamiltonian cannot be used as A u s /a is not a small quantity compared with 
the energy cutoff u s /a. To describe the low energy physics at A = 7r, it is necessary to choose a gauge with the vector 
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potential along the rungs of the ladder, so that the interchain hopping reads: 

Hhop. = ^ 

3 

Applying bosonization to (IS32I) . we obtain the following representation for interchain hopping: 


Hhop. — 


n 

2ira 


dx cos V2(j> c L-<v^(e.+0.) + e -iV2(e s -4>s) + e iV 2 (e s +<t> a ) + e 


which can be rewritten in terms of SU(2) 1 Wess-Zumino-Novikov-Witten (WZNW) currents[S19J: 

Hhop. = Q j dx cos V2(p c (J R + J V L ). 

The resulting Hamiltonian H c + H s + Hh op . can be treated in mean-field theory S20l - S23 |. giving: 

TT _ TtMF I ttMF 

MF ■^c t i 


H. 


MF 


H. 


MF 


f dx 

J 2tt 

2nu s 


u c K c (ttU c ) 2 + -^(d x (j ) c ) 2 


+ — dx cos V2d> Ci 
na J 

J dx(J R -Jr + Jl- j l) + h s J dx(J y R + J y L ), 


(532) 

(533) 

(534) 

(535) 

(536) 

(537) 


where: 


— = 8 n(J y R + J y L ) s ,MF, 

na 

h s — 8f2(cos V 24 >c) c ,mf- 


(S38) 


Using a rotation around the x axis, = J*, J* = — and applying abelian bosonization [S24| , we rewrite: 


ttMF _ / 

Hs ~ / 2 ^ s 


dx 


(7rfl s ) 2 + ( d(j) s f 


h s 


tV2 . 


d x (j) s dx, 


which allows us to write: 


U<w= E <■£> = (J V R +jd = - 


tV? 


v—R,L 


2ttu s ’ 


(S39) 


(S40) 


and allows us to solve (1S38I) with h s ~ H 2 and g c ~ fi 3 . We obtain a gap in the total density excitations, A c ~ fl 2 , while 
the antisymmetric modes remain gapless and develop an incommensuration. To characterize the incommensuration, 
we need to detail the rotation of the SU(2), WZNW currents and primary fields. After shifting <j) s — > <j) s + h ‘% , we 
find: 


1 


rV2 


d x <fi s =-E 

xYs 2na ^ 


irV2(9,+r'<p s )+irr' 


r,r / =± 


E sin = e sin f 9 S 


r=± 


r=± 


&x4*s 


and: 


E cos >/ 2 ( 0 s + =-^ 7 = 

r=± ^2 


sin \PlQ s = sin -\/20s 

cos \[2Q S = sin 

sin = — cos V20 s 

cos \Pl(j) s = cos ( V2(j) s + ] 

V Us J 


2ttu s ’ 


(541) 

(542) 

(543) 

(544) 

(545) 

(546) 

(547) 
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Since we have: 


3±ti) 

a*(x) 

0 C BOW 

we find, after the rotation: 


— ^ sin y/2(0 S + rcj) s ) + 
na z —f 

r=± 


2 n(~y 

TTCL 


sin V29 S , 


1 (-1V 

dx4> s + 


nV2 

(-l)J r 
-(cos V 2^ c ) cos 


(cos V2<j) c ) sin V2cj) s , 
V2tj> s 


j±(j) = — V' sin y/2 f8 S + r<j) s + r-^-\ + ——— sin y/28 s , 
na ■*—' V u s J na 

r=± x x 

= — V cos\/2 [ 0 S + r(j) s + r—— ] — -—— (cos V2(j> c ) cos\/2 6 S , 
7Ta f V U s J 7TCL 

r—zt x x 

°scw = ^ — (cos\/2^ c )cos (V 24 > s + 

TTCL \ 'Ug y 


(548) 

(549) 

(550) 

(551) 

(552) 

(553) 


so that: 


{j±{j)j±(j')) 

<* z (jV(/)> 

{' 0 C BOw(j)0 C BOw{j ')) 


1 


2 71-2 (j ~ j') 2 
1 

2 t T 2 {j-j') 2 


COS 


cos 


hs{j ~f) 
Us 


-p-p-cos 

b — j I 


U - j'1 

1/ -j'l 


(S54) 


(S55) 


(S56) 


We see that an incommensuration of wavevector p(fi) = h s /u s develops in the k ~ 0 component of the rung current 
and density wave correlations. Since p(fi) ~ id 2 (see Fig. lS2l) . the incommensuration increases with interchain hopping. 
The Fourier transform of the k ~ 0 component behaves as |k — p(Sl)\ + \k + p(f2)|, i. e. it is constant for |fc| < p(f2) 
and linear in k for |fc| > p(fl). 



FIG. S2: A graph of the incommensuration p(fl) = h s /u s at A = 7r obtained from numerical data (L=64 in PBC). The blue 
dots are the value of S(k = 0), while the black dots correspond to the slope discontinuity k = p(fi). The red and green lines 
are quadratic fits. In the inset we show S(k = 0), open blue dots, together with the Ksp(Q)/n, black soid dots, with the choice 
K* = 1 

Since the boson annihilation operators do not correspond to primary fields of the SU( 2)i WZNW model, we cannot 
directly derive their expression from SU( 2) symmetry. However, since e l9 “l^ 2 has conformal dimensions (1/16,1/16) 
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its expression in terms of 6 S and <j> a has to be a sum of operators of conformal dimensions (1/16,1/16). A general 
expression is: 


eV 2 = Aoe-ss + + A 2 e l ^+A 3 e 


Moreover, the operator product expansion: 


.w 

e' va 


75 


X — x' 
a 


1/4 




(S57) 


(S58) 


has to be satisfied, so, for A = n, w e must have {A*,, Aj} = 6kje lk 5 which can be satisfied by writing Ak as the 
product of a 4 x 4 Dirac matrix [S25j by a phase e~ lk T. In the case of A = 7r we obtain the correlator: 


(e 75 


3s(») 4 Qs(x') 

e 75 


) = 


(. “ ,.y 

2 + 2 cos —(x — x') 

\\x-x'\J 

u s 


(S59) 


where the cos results from the correlation of the This implies that the momentum distribution is: 


where: 


r(*) = 


Jdxe^(e q ’ < % { - x) e l< % { - 0) ){e 6B ^ ) e ) 
= I dxe i( ' k+7Ta ' >x (e _i (x) e* (0) ) ( — 


2 + 2 cos 


h s x 


1 / h s \ 1 / h s 

= vik + 7rer) + —i>[ k + na -\ -+ —v\ k + na - 

2 V Us J 2 V Uc 



dxe ikx {e~ i ^ (x) e i ^ 0) ) 



(560) 

(561) 

(562) 


(S63) 


Since the width at half maximum of the Lorentzian-shaped graph of the function v scales as A c /u c ~ D 2 and 
h s /u s ~ D 2 , the graph of n(k ) can either comprise 3 peaks or a single broad peak centered in 7 ra depending on the 
dimensionless ratio h s u c /(u s A c ). 

When A < 7T, we choose a gauge such that: 

H = -*£( 6 L ei(A ~* )ff w+H.c.) 

3,0- 

+D^(-l)^t CT 6 j) _ CT , (S64) 

3,° 

and we define <5A = A — 7r. The mean field Hamiltonian becomes Hmf = H^ F + H^ IF with: 

Hf F = ^ J dx{ Jr ■ Jr + J L • j l) +h s J dx{J v R + f}) + ^ J dx(J z R - J£), (S65) 


and H^ f unchanged. We now have to make a different rotation around x for the right moving and the left moving 
current: 


// 


MF 


J V H 

= sin tpJft + cos tpJ R 

(S66) 

Jr 

= — cos ipJ R + sin +J R 

(S67) 

J V l 

= — sin ipj^ + cos (pJ[ 

(S68) 

T z 

J L 

= — cos </> J| H— sin <pJ R 

(S69) 

dx 
——u 
27r 

S (7rfl s) 2 + (d(j) s ) 2 - J d x <p s dx, 

(S70) 


To find; 

















where h s ( A) = y/h% + u 2 s (5X/a) 2 . We still have (+ J|) = — so the mean-field equations remain the same. We 
also find: 


sin V20 S = 


h s . r-~ udX/a 

sin v 20 s + - — 77 T- cos 


M A) 

( 


M A) 


(V^ s + ^x 


cos-\/20 s = sin ( V2(f> s + ——-x 


sin \J~2(f> s = — cos V20 S 


cos V2 (j) s = 


M A) 


cos 


y/2<j> s 


h s ( A) 


The staggered part of the rung current correlations becomes: 


(j±{x)j±(0)) 


uSX/a 

MA) 


(-)*/“ hl+(^y C os^x 


sin 


V29 s 


(571) 

(572) 

(573) 

(574) 

(575) 


M A) 2 

so that the Fourier transform will present peaks at k = ^ and k = — ± h s {X)/u s . When <5 A is increased, the two peaks 
at k = — ± become dominant, and we crossover to the behavior already discussed for weak A. In the case of 

a u s 1 J 

S(k), the peak at k = n is not split as A is reduced. If we look at the BOW c correlations, a peak at k = n appears, 
and becomes the dominant peak when <5A is increased. Using the rotation (IS66I) we can also obtain the antisymmetric 
density correlations as: 

f u s SX \^ 

K a ) 1 


^{d x <t> s {x)dM o)} = ^2 


hi 


■ cos( h s (X)/u s x) 


(S76) 


h2 sW h* + (2*py 

When <5A increases, this expression crosses over to the 1/(27r 2 x 2 ) which was obtained at small A. S(k) now presents 
a change of slope at |fc| = h s (X)/u s . 

As to the k ~ 0 component of the rung current, since it is proportional to it becomes Jr + under the 

rotation, and the correlator becomes: 

'MA),. 


1 

2?r 2 {j ~j') 2 


cos 


-(3-3')^ 


(S77) 
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